Three common mistakes in statistics and how to avoid them

class: center, middle, inverse, title-slide

.title[
# Three common mistakes in statistics and how to avoid them
]
.author[
### Elizabeth Pankratz
]
.institute[
### Department of Psychology<br>The University of Edinburgh
]

---

class: middle

## But first:

[mentimeter link]

[mentimeter qr code]

---

.pull-left[
## The mistake
]

.pull-right[
## How you'll avoid it
]

.pull-left[
<svg aria-hidden="true" role="img" viewBox="0 0 640 512" style="height:1em;width:1.25em;vertical-align:-0.125em;margin-left:auto;margin-right:0.2em;font-size:inherit;fill:currentColor;overflow:visible;position:relative;"><path d="M392.8 1.2c-17-4.9-34.7 5-39.6 22l-128 448c-4.9 17 5 34.7 22 39.6s34.7-5 39.6-22l128-448c4.9-17-5-34.7-22-39.6zm80.6 120.1c-12.5 12.5-12.5 32.8 0 45.3L562.7 256l-89.4 89.4c-12.5 12.5-12.5 32.8 0 45.3s32.8 12.5 45.3 0l112-112c12.5-12.5 12.5-32.8 0-45.3l-112-112c-12.5-12.5-32.8-12.5-45.3 0zm-306.7 0c-12.5-12.5-32.8-12.5-45.3 0l-112 112c-12.5 12.5-12.5 32.8 0 45.3l112 112c12.5 12.5 32.8 12.5 45.3 0s12.5-32.8 0-45.3L77.3 256l89.4-89.4c12.5-12.5 12.5-32.8 0-45.3z"/></svg>
**A common R programming mistake:**
Letting R treat all variables that consist of numbers as numeric.
]

.pull-right[
]

.pull-left[
<svg aria-hidden="true" role="img" viewBox="0 0 448 512" style="height:1em;width:0.88em;vertical-align:-0.125em;margin-left:auto;margin-right:0.2em;font-size:inherit;fill:currentColor;overflow:visible;position:relative;"><path d="M160 80c0-26.5 21.5-48 48-48h32c26.5 0 48 21.5 48 48V432c0 26.5-21.5 48-48 48H208c-26.5 0-48-21.5-48-48V80zM0 272c0-26.5 21.5-48 48-48H80c26.5 0 48 21.5 48 48V432c0 26.5-21.5 48-48 48H48c-26.5 0-48-21.5-48-48V272zM368 96h32c26.5 0 48 21.5 48 48V432c0 26.5-21.5 48-48 48H368c-26.5 0-48-21.5-48-48V144c0-26.5 21.5-48 48-48z"/></svg> 
**An advanced statistical mistake:**
Modelling categorical, ordinal data as if it were numeric.
]

.pull-right[
]

.pull-left[
<svg aria-hidden="true" role="img" viewBox="0 0 384 512" style="height:1em;width:0.75em;vertical-align:-0.125em;margin-left:auto;margin-right:0.2em;font-size:inherit;fill:currentColor;overflow:visible;position:relative;"><path d="M192 32c17.7 0 32 14.3 32 32V199.5l111.5-66.9c15.2-9.1 34.8-4.2 43.9 11s4.2 34.8-11 43.9L254.2 256l114.3 68.6c15.2 9.1 20.1 28.7 11 43.9s-28.7 20.1-43.9 11L224 312.5V448c0 17.7-14.3 32-32 32s-32-14.3-32-32V312.5L48.5 379.4c-15.2 9.1-34.8 4.2-43.9-11s-4.2-34.8 11-43.9L129.8 256 15.5 187.4c-15.2-9.1-20.1-28.7-11-43.9s28.7-20.1 43.9-11L160 199.5V64c0-17.7 14.3-32 32-32z"/></svg>
**A foundational statistical mistake:**
Interpreting a significant *p*-value as proof that an effect exists.
]

.pull-right[
]

---
class: inverse, middle, center

# The data we'll use

---

## The SMARVUS dataset (Terry et al., 2023)

.center[SMARVUS = **S**tatistics and **M**athematics **A**nxieties and **R**elated **V**ariables in **U**niversity **S**tudents]

.pull-left[
A survey of *n* = 18,841 students (mostly Psychology UGs) from 35 countries.

Students rated their anxiety from 1 (no anxiety) to 5 (a great deal of anxiety) in scenarios like:

- Studying for a statistics test.

- Interpreting the meaning of a table in a journal article.

- **Going to ask my statistics teacher for individual help with material I am having difficulty understanding.**
]

.pull-right[

]

---

## Why Likert scale ratings are not numeric

.center[
![:scale 70%](data:image/png;base64,#imgs/rating-scale.png)
]

---

## And yet...

.center[
![:scale 50%](data:image/png;base64,#imgs/ReederEtAl2017-JML.png)
]

Figure 2 from Reeder et al. (2017), published in Journal of Memory and Language!

---

## By default, R will keep numeric-looking variables numeric

If we allow R's default behaviour, then we can do naughty things with categorical variables:

``` r
mean(anx$rating)
```

```
## [1] 2.868054
```

But if we store these variables as factors, the naughty things become impossible (yay!):

``` r
anx <- anx |>
  mutate(rating = factor(rating))

mean(anx$rating)
```

```
## [1] NA
```

---

# LOs

---
class: inverse, middle, center

# Modelling an ordinal variable

### The .mono-white[polr()] express

---

## Fit an ordinal regression model

``` r
library(MASS)        # MASS contains the polr() function
                     # (polr = Proportional Odds Logistic Regression)

anx_fit1 <- polr(
  rating ~ 1,        # intercept-only model, to start
  data = anx, 
  Hess = TRUE,       # required if we want to use summary()
  method = 'probit'  # more on this in a moment
)
```

---

## Fit an ordinal regression model

``` r
summary(anx_fit1)
```

```
## Call:
## polr(formula = rating ~ 1, data = anx, Hess = TRUE, method = "probit")
## 
## No coefficients
## 
## Intercepts:
##     Value    Std. Error t value 
## 1|2  -0.8420   0.0157   -53.7268
## 2|3  -0.1678   0.0138   -12.1462
## 3|4   0.3833   0.0141    27.1512
## 4|5   1.0339   0.0168    61.6193
## 
## Residual Deviance: 26596.28 
## AIC: 26604.28
```

---

## What do those `Intercepts` mean?

???

- imagine that there's some underlying continuous normal distribution of anxiety, assumed standard normal [show normal distrib]
- ppl with high anxiety are more likely to give high responses, ppl with low anxiety more likely to give low responses (could do emojis relating to anxiety: <svg aria-hidden="true" role="img" viewBox="0 0 512 512" style="height:1em;width:1em;vertical-align:-0.125em;margin-left:auto;margin-right:auto;font-size:inherit;fill:currentColor;overflow:visible;position:relative;"><path d="M464 256A208 208 0 1 0 48 256a208 208 0 1 0 416 0zM0 256a256 256 0 1 1 512 0A256 256 0 1 1 0 256zm177.6 62.1C192.8 334.5 218.8 352 256 352s63.2-17.5 78.4-33.9c9-9.7 24.2-10.4 33.9-1.4s10.4 24.2 1.4 33.9c-22 23.8-60 49.4-113.6 49.4s-91.7-25.5-113.6-49.4c-9-9.7-8.4-24.9 1.4-33.9s24.9-8.4 33.9 1.4zm40-89.3l0 0 0 0-.2-.2c-.2-.2-.4-.5-.7-.9c-.6-.8-1.6-2-2.8-3.4c-2.5-2.8-6-6.6-10.2-10.3c-8.8-7.8-18.8-14-27.7-14s-18.9 6.2-27.7 14c-4.2 3.7-7.7 7.5-10.2 10.3c-1.2 1.4-2.2 2.6-2.8 3.4c-.3 .4-.6 .7-.7 .9l-.2 .2 0 0 0 0 0 0c-2.1 2.8-5.7 3.9-8.9 2.8s-5.5-4.1-5.5-7.6c0-17.9 6.7-35.6 16.6-48.8c9.8-13 23.9-23.2 39.4-23.2s29.6 10.2 39.4 23.2c9.9 13.2 16.6 30.9 16.6 48.8c0 3.4-2.2 6.5-5.5 7.6s-6.9 0-8.9-2.8l0 0 0 0zm160 0l0 0-.2-.2c-.2-.2-.4-.5-.7-.9c-.6-.8-1.6-2-2.8-3.4c-2.5-2.8-6-6.6-10.2-10.3c-8.8-7.8-18.8-14-27.7-14s-18.9 6.2-27.7 14c-4.2 3.7-7.7 7.5-10.2 10.3c-1.2 1.4-2.2 2.6-2.8 3.4c-.3 .4-.6 .7-.7 .9l-.2 .2 0 0 0 0 0 0c-2.1 2.8-5.7 3.9-8.9 2.8s-5.5-4.1-5.5-7.6c0-17.9 6.7-35.6 16.6-48.8c9.8-13 23.9-23.2 39.4-23.2s29.6 10.2 39.4 23.2c9.9 13.2 16.6 30.9 16.6 48.8c0 3.4-2.2 6.5-5.5 7.6s-6.9 0-8.9-2.8l0 0 0 0 0 0z"/></svg>, <svg aria-hidden="true" role="img" viewBox="0 0 512 512" style="height:1em;width:1em;vertical-align:-0.125em;margin-left:auto;margin-right:auto;font-size:inherit;fill:currentColor;overflow:visible;position:relative;"><path d="M256 48a208 208 0 1 0 0 416 208 208 0 1 0 0-416zM512 256A256 256 0 1 1 0 256a256 256 0 1 1 512 0zM168 320c-13.3 0-24 10.7-24 24s10.7 24 24 24h8V320h-8zm40 48h32V320H208v48zm96 0V320H272v48h32zm32 0h8c13.3 0 24-10.7 24-24s-10.7-24-24-24h-8v48zM168 288H344c30.9 0 56 25.1 56 56s-25.1 56-56 56H168c-30.9 0-56-25.1-56-56s25.1-56 56-56zm-23.6-80a32 32 0 1 1 64 0 32 32 0 1 1 -64 0zm192-32a32 32 0 1 1 0 64 32 32 0 1 1 0-64z"/></svg>)
- so to estimate how different anxiety levels translate to different responses on the 1--5 scale, we draw thresholds on that distribution [add thresholds]
- ppl with anxiety in this bin will respond with 1, in this bin with 2, etc.
- and those thresholds, the cutpoints btwn ratings, are the intercepts.
- [show intercept estimates, put thoes same numbers on the thresholds]
- normal distribution assumption is from method = probit. other methods assume other underlying distributions, but the idea of thresholds is the same.

---
class: inverse, middle, center

# "Preregistration" activity

### The effect of gender on stats anxiety

---

#### How will a student's gender affect ratings for "Going to ask my statistics teacher for individual help with material I am having difficulty understanding"?

---

Activity is part of the handout.

???

The puzzle pieces:

- The threshold values will stay in the same place.
- If you move the normal distribution left, then more probability density will end up toward the lower end of the scale (less anxiety). If you move the distribution right, more probability density will be toward the higher end of the scale (more anxiety).
- Each gender gets one probability distribution over anxiety values.
- Here's the one for Female/Woman.
- I want you to "preregister" a guess, before we look at the outcome of the model, whether the distribution for the Male/Men group will be shifted left or right from the Female one, and whether the distribution for the Another Gender group will be shifted left or right from the female one. Based on this data.
- Draw it on the plot on the handout, if you want.

---

## How does gender affect anxiety ratings?

.center[
<img src="data:image/png;base64,#demo_files/figure-html/all-gender-normals-1.png" width="864" style="display: block; margin: auto;" />
]

---

``` r
anx_fit2 <- polr(
  rating ~ gender,
  data = anx,
  method = 'probit',
  Hess = TRUE
)
summary(anx_fit2)
```

```
## Call:
## polr(formula = rating ~ gender, data = anx, Hess = TRUE, method = "probit")
## 
## Coefficients:
##                        Value Std. Error t value
## genderMale/Man       -0.3280    0.03015 -10.880
## genderAnother Gender  0.4846    0.11992   4.041
## 
## Intercepts:
##     Value    Std. Error t value 
## 1|2  -0.9045   0.0169   -53.5402
## 2|3  -0.2246   0.0150   -14.9847
## 3|4   0.3318   0.0151    21.9158
## 4|5   0.9889   0.0176    56.2958
## 
## Residual Deviance: 26456.99 
## AIC: 26468.99
```

---

# LOs

---
class: inverse, middle, center

# Interpreting *p*-values

---

## Are the effects of `gender` significant?

```
Coefficients:
                       Value Std. Error t value
genderMale/Man       -0.3280    0.03015 -10.880
genderAnother Gender  0.4846    0.11992   4.041
```

No *p*-values in the model summary.
But it's common practice to compare `t value` to a standard normal distribution (treating them like *z* scores).

``` r
pnorm(abs(-10.880), lower.tail = FALSE) * 2
```

```
## [1] 1.43563e-27
```

``` r
pnorm(abs(  4.041), lower.tail = FALSE) * 2
```

```
## [1] 5.322376e-05
```

Since both *p*-values are below 0.05:

- we CAN reject the null hypothesis that gender has no effect on ratings.
- we CANNOT conclude that there really is an effect of gender.

---

#### Why doesn't a significant *p*-value mean that there really is an effect?

Because we can also get significant *p*-values when there really is *no* effect.

.pull-left[
No difference in reality:

<img src="data:image/png;base64,#demo_files/figure-html/true-skew-probdist-1.png" width="504" style="display: block; margin: auto;" />
]

.pull-right[
A possible random sample (*n* = 50 per group):

]

---

``` r
sim_fit <- polr(rating ~ group, data = simdat, method = 'probit', Hess = TRUE)
summary(sim_fit)  
```

```
## Call:
## polr(formula = rating ~ group, data = simdat, Hess = TRUE, method = "probit")
## 
## Coefficients:
##                Value Std. Error t value
## groupGroup B -0.4479     0.2229  -2.009
## 
## Intercepts:
##     Value   Std. Error t value
## 1|2 -0.3244  0.1677    -1.9344
## 2|3  0.1738  0.1664     1.0447
## 3|4  0.9263  0.1879     4.9283
## 4|5  1.4570  0.2333     6.2444
## 
## Residual Deviance: 267.4865 
## AIC: 277.4865
```

``` r
pnorm(abs(-2.009), lower.tail = FALSE) * 2
```

```
## [1] 0.04453713
```

---

.pull-left[
## The mistake
]
.pull-right[
## How you'll avoid it
]

.pull-left[
<svg aria-hidden="true" role="img" viewBox="0 0 640 512" style="height:1em;width:1.25em;vertical-align:-0.125em;margin-left:auto;margin-right:0.2em;font-size:inherit;fill:currentColor;overflow:visible;position:relative;"><path d="M392.8 1.2c-17-4.9-34.7 5-39.6 22l-128 448c-4.9 17 5 34.7 22 39.6s34.7-5 39.6-22l128-448c4.9-17-5-34.7-22-39.6zm80.6 120.1c-12.5 12.5-12.5 32.8 0 45.3L562.7 256l-89.4 89.4c-12.5 12.5-12.5 32.8 0 45.3s32.8 12.5 45.3 0l112-112c12.5-12.5 12.5-32.8 0-45.3l-112-112c-12.5-12.5-32.8-12.5-45.3 0zm-306.7 0c-12.5-12.5-32.8-12.5-45.3 0l-112 112c-12.5 12.5-12.5 32.8 0 45.3l112 112c12.5 12.5 32.8 12.5 45.3 0s12.5-32.8 0-45.3L77.3 256l89.4-89.4c12.5-12.5 12.5-32.8 0-45.3z"/></svg>
**A common R programming mistake:**
Letting R treat all variables that consist of numbers as numeric.
]
.pull-right[
When you know a variable is categorical, tell R that using `factor()`.
]

.pull-left[
<svg aria-hidden="true" role="img" viewBox="0 0 448 512" style="height:1em;width:0.88em;vertical-align:-0.125em;margin-left:auto;margin-right:0.2em;font-size:inherit;fill:currentColor;overflow:visible;position:relative;"><path d="M160 80c0-26.5 21.5-48 48-48h32c26.5 0 48 21.5 48 48V432c0 26.5-21.5 48-48 48H208c-26.5 0-48-21.5-48-48V80zM0 272c0-26.5 21.5-48 48-48H80c26.5 0 48 21.5 48 48V432c0 26.5-21.5 48-48 48H48c-26.5 0-48-21.5-48-48V272zM368 96h32c26.5 0 48 21.5 48 48V432c0 26.5-21.5 48-48 48H368c-26.5 0-48-21.5-48-48V144c0-26.5 21.5-48 48-48z"/></svg> 
**An advanced statistical mistake:**
Modelling categorical, ordinal data as if it were numeric.
]
.pull-right[
Apply and interpret ordinal regression models (e.g., `polr()` from `MASS`).
]

.pull-left[
<svg aria-hidden="true" role="img" viewBox="0 0 384 512" style="height:1em;width:0.75em;vertical-align:-0.125em;margin-left:auto;margin-right:0.2em;font-size:inherit;fill:currentColor;overflow:visible;position:relative;"><path d="M192 32c17.7 0 32 14.3 32 32V199.5l111.5-66.9c15.2-9.1 34.8-4.2 43.9 11s4.2 34.8-11 43.9L254.2 256l114.3 68.6c15.2 9.1 20.1 28.7 11 43.9s-28.7 20.1-43.9 11L224 312.5V448c0 17.7-14.3 32-32 32s-32-14.3-32-32V312.5L48.5 379.4c-15.2 9.1-34.8 4.2-43.9-11s-4.2-34.8 11-43.9L129.8 256 15.5 187.4c-15.2-9.1-20.1-28.7-11-43.9s28.7-20.1 43.9-11L160 199.5V64c0-17.7 14.3-32 32-32z"/></svg>
**A foundational statistical mistake:**
Interpreting a significant *p*-value as proof that an effect exists.
]
.pull-right[
Understand that significant *p*-values can arise even if no effect exists.
]

---

## TODO References

Reeder et al. 2017

Terry et al. 2023